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Abstract 

Coupling of the Reynolds- averaged Navier-Stokes equations, rigid-body dynamics, and a pitch at- 
titude control law is demonstrated in two- and three-dimensions. The application problem was the 
separation of a canard-controlled store from an open-flow rectangular cavity bay at a freestream Mach 
number of 1.2. The transient flowfield was computed using a diagonal scheme in an overset mesh frame- 
work, with the resultant aerodynamic loads used as the forcing functions in the nonlinear dynamics 
equations. The proportional and rate gyro sensitivities were computed apriori using pole placement 
techniques for the linearized dynamical equations. These fixed gain values were used in the controller 
for the nonlinear simulation. Reasonable comparison between the full and linearized equations for a 
perturbed two-dimensional missile was found. Also in two-dimensions, a controlled store was found 
to possess improved separation characteristics over a canard-fixed store. In three-dimensions, trajec- 
tory comparisons with wind-tunnel data for the canard-fixed case will be made. In addition, it will 
be determined if a canard-controlled store is an effective means of improving cavity store separation 
characteristics. 


Nomenclature 


C m 

coefficient of pitching moment 

DOF 

degree of freedom 

I 

body inertia 

Ka 

proportional gain 

K r 

rate gyro sensitivity 

L 

characteristic length 

in 

body mass or stage number 

M 

Mach number 

Q 

dynamic pressure 

Q 

vector of dependent variables 

Re 

Reynolds number 

s 

Laplace operator 
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f time 

T absolute' temperature 

u. ik w Cartesian velocity components 

x, y. z Cartesian physical space 

a angle of attack 

effector deflection angle 
p density 

rr shear layer spreading rate parameter 

r canard servo characteristic time 

0 body attitude 

oj vorticity 

( ) time derivative 

Superscripts 

n time step 

Subscripts 

c commanded, canard 

t tail 

oc freestream quantity 


Introduction 

The design of high-performance aircraft has in the past typically been compartmentalized by dis- 
cipline. most commonly structures, fluids, and controls. However, as performance demands escalate, 
aircraft systems have become increasingly interrelated. 1 ' 2 Therefore, there is a need to investigate the 
optimization of the aircraft in its entirety, not simply by evaluation of its sub-systems. 

This type of multidisciplinary analysis is currently accomplished with simplified physical models, 
such as panel flow and modal structural codes. However, in critical regions of the flight envelope those 
linear methods can fail, leading to the necessity for higher-order models.' 1 ' 1 Obviously, this increased 
physical fidelity comes with a high computational price, and hence fewer design cycles are permitted 
as compared to the linear methods. 

Providing a capability for the simulation of the nonlinear interaction of fluids and rigid body motion 
will be useful in several ways. First, the simulation could be used to validate the implementation of a 
control law derived using conventionally determined sensitivity coefficients. Although one must be wary 
of results obtained numerically, simplifications used to compute these force and moment derivatives 
will not be present in a Navier-Stokes simulation. Second, the coupled simulations could be used to 
develop a control law where nonlinear effects are important, using the computed aerodynamic forces and 
moments instead of tabulated empirical relationships. 0 Hence, computational design and prototyping 
of the aircraft control system offers the capability of reducing aircraft design cycle cost and enhancing 
safety as well as improving performance. 

The effort, documented here begins to address the interaction of the disciplines of fluids, rigid 
body dynamics, and controls in nonlinear flight regimes. In order to assess the accuracy of these 
initial computations, a problem which could be compared against analytic and experimental results 
was chosen: the cavity store separation problem. Dix and Dobson's 6 recent experimental study of 
the separation of stores from cavity bays will be used as a basis for comparison. These wind-tunnel 
tests determined the trajectory of an uncontrolled missile from a rectangular cavity and will be used 
to validate the canard-fixed computation. 

Previous computational efforts have shown that the component problems of cavity flows' ' 9 and 

uncontrolled store separation 10, 11 can be solved with reasonable accuracy using overset grid methods. 
It should be noted that inviscid solutions to the cavity store separation problem, with bodies in relative 
motion, have been obtained on tetrahedral meshes. 12 However, although these unstructured methods are 
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geometrically powerful, a means of consistently accounting for viscous effects is currently lacking. Here, 
r lit' combined problem of viscous cavity How. rigid body dynamics, and automatic control technique* 
is addressed in a general manner using the overset mesh framework. 

The following sections discuss the approach used to solve the coupled system and the results ob- 
tained for several two- and three-dimensional cases. Comparisons of numerical results are math' against 
linearized or experimental results where available. 

Approach 

The procedure used to solve these coupled problems is shown schematically in Fig. 1. The process 



Fig. 1: Overall coupled system approach 

begins with the generation of the blocked mesh system while holding the relative position of the 
geometry in the initial position. After establishing the initial grid connectivity information, the solution 
of the flowfield can begin. Obviously, depending on the problem being addressed, this fixed grid solution 
may be steady or possess a bounded envelope. When this fixed-geometry flow solution is satisfactory, 
the motion of the body commences. 

The integrated aerodynamic loads acting on the body are provided to the 6-DOF portion of the 
code along with the body state vector. Using this information, the body position and attitude is then 
integrated one time step, and kinematic constraints are applied. The controller also uses the body state 
to compute new effector settings, which are then applied to all relevant grids. Finally, since the grids 
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have changed relative positions, the intergrid communication is re-established. This process repeats 
for oat'll time stop until tlio simulation is romploto. A dot.ailod description of tin 1 component processes 
is given below. 


Flow Solver 

The fluid field was computed via the Reynolds-averaged Navier-Stokes (RANS) equations using the 
diagonal scheme of Pulliam and Chaussee 14 implemented in the overset grid framework of Benek. Bull- 
ing, and Steger. 14 The equations were integrated through Euler implicit time marching and second-order 
spatial differencing with viscous wall conditions specified as no-slip, zero normal pressure gradient, and 
adiabatic. Information transfer across overset mesh boundaries was implemented using trilinear inter- 
polation of the dependent variable vector, Q = [ p , pu, pv, pw, e ] T . The flow solver cost is 13/ns/cell/step 
on a Cray Y-MP. 


Turbulence Model 


In these computations, the slowly time-varying component of the flow is resolved, while rapid fluc- 
tuations are modelled. The algebraic turbulence model of Baldwin and Lomax, 15 as implemented by 
Buning, 16 is described below using a flow in the (x,y) plane. 

The wall-bounded flows use the original Baldwin- Lomax model, with the addition of a variable 
Fmax cutoff. The grids are chosen such that a unique wall distance is readily available. 

The cavity aperture spanning shear layer uses an eddy viscosity developed using F(y) = y |w|, as 
suggested by Baldwin and Lomax 10 for wake regions. This results in 


y wake 


OC 


cu- 


Vmar u dif 


c wk m 

\ jkq max ) 



where specification of C wk is discussed below and the velocity difference is modified to be half the total 
velocity difference between the streams in the specified shear layer region 

u d ,f = }/{u 2 + v 2 ) max - >/(« 2 + « 2 ) Mm „ 


The free shear layer model is now given by 


pi — pFCcpCwk 


u if 


Mr 


(1) 


after dropping the Klebanoff intermittency function. Unmodified model constants K = 0.0168 and 
C C p = 1.6 are used. 

The magnitude of the eddy viscosity in the free shear layer model can be altered with C wk . Esti- 
mation of the proper value of C wk begins by using Gdrtler’s shear layer solution: 
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where u\ and 112 are the velocities of the slow and fast streams and £ is the similarity coordinate. The 
spreading parameter a is inversely related to the spreading rate, db/dx , where b is a measure of the 
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shrar layrr width, ['hr value of tin 1 spreading parameter when the velocity of one of the stream^ i- 
zer< > is <i\ |. 

Gdrtler's solution can he used to determine the maximum vorticity magnitude as follows: 
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Now. using PrandtFs mixing length assumption and scaling laws for jet boundaries, eddy viscosity 
can also be expressed as 

/if OC pi 2 \u\ 


oc pi 2 


> Au 


= KopbAu 


( 2 ) 


where Kq = 


T 



Setting Eqs. 1 and 2 equal results in 


G — ( &oKCs C p \/^0 

and only ctq remains to be specified. Empirical estimates of oq range from 9.0 to 13. 5. 9 For this series 
of cavity flow efforts (tq was set to 11.0, resulting in a value of C w ^ = 1.91. 


Grid Generation and Communication 

Computation of the loads generated by cavity flows requires accurate representation of the geometry 
as well as the flowfield. Typically, a significant effort in grid generation is required before flow analysis 
can begin. Since matching zone faces are not required for the overset method used here, recent advances 
in algebraic 1 ' and hyperbolic 18 methods can be used. Hyperbolic grid generation, which gives good 
spacing and orthogonality control, was used for the wall-bounded regions, while algebraic grids were 
used in shear flow regions. Advantages of this type of grid system include straightforward specification 
of the turbulent regions and allowance for independent refinement of each zone. This topology also 
permits the re-use of meshes for configuration studies. 

Exchange of flow information is accomplished using a domain connectivity function, with the donor- 
receiver relationship established at each time step using an efficient technique. 19 Although the cost of 
re-establishing intergrid communication is problem dependent, the computational expense is generally 
half of that used by the diagonalized flow solver. The initial location for the hole-cutters was specified 
using a graphical interface, 20 after which the movement of the grids and hole cutters were updated 
automatically. An example of the overset mesh topology used is shown in Fig. 2. which shows the 
configuration at an instant in the separation process. 


Kinematics and Rigid Body Dynamics 

Although details of the six-DOF dynamical equations can be found elsewhere 21 briefly the rotational 
dynamics is described by Euler's equations of motion which aligns the xyz coordinates with the body 
principal axes at the center of gravity. For instance, for a rotation 6 about an axis A, Euler parameters 
can be specified as 

I, . 0 . e . 9 ep 

€= |Aisin-, A2S111-, A3 sin-, cos 
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Fig. 2: Coarsened grids after release 


These Euler parameters, integrated according to the rotational body dynamics, are updated and stored 
for each grid. 10, 19 Kinematic constraints can be imposed during the ejection process or for restricted- 
DOF simulations. In addition, the assumption of rigid-body dynamics eliminates the need to store the 
component grids for all time, since the Euler parameters may be used to compute grid attitude from 
the initial position. 

Rotating effectors were implemented by summing the relative commanded effector and body angular 
velocities, u). Integration of u gives the proper effector attitude relative to the initial grid positions. 
The hinge line location is updated according to the attitude of the body. Storage of the hinge line and 
Euler parameters associated with each grid allows nesting of parent-child bodies to an arbitrary level 
without modification to the grid communication and support software. 

Pitch Attitude Control Law 

The fourth-order system shown in Fig. 3 was used as the system model for both the two- and 



Fig. 3: 2-D missile: Block diagram 


three-dimensional missile cases. In Fig. 3 the plant and servo can be expanded as 


e_ 
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ds ~F e 
as 2 + bs + c 


;/ = 


l 

r 


where the servo time constant was taken as r = J? The plant coefficients were determined from 
the governing equation; they are composed of geometric and flow information, as well as the stability 
derivatives. The stability derivatives were determined from linearized supersonic airfoil theory in the 
two-dimensional cases and from direct computation in three-dimensions. Pole placement techniques 
and linearized system time response were used to determine the proportional. K a , and the rate gyro, 
K r „ sensitivities. 
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Results and Discussion 


The method described above was applied ill two-dimensions to a perturbed one degree-of-freedom 
ease and a tliree degree-of-freedom eavity store separation process. The three-dimensional results are 
in progress. 


One Degree-of- Freedom Simulation 

In order to provide some measure of validation, a two-dimensional, 1-DOF simulation was imple- 
mented at a flight altitide of 45, 000 ft. Comparison of the linearized and the coupled nonlinear system 
response for a small perturbation allowed assessment of the basic methodology. 


Linearized Dynamics 

The equation of motion, M y = I yy 9, after neglecting the pitch-damping term, is 

My = qo e {[(C, a Sd) c - (Ci a Sd)t]9 + ( C, a Sdfi) c } 

where the body was pinned at the location of the center of gravity for this analysis. 

Expressing the system in a state variable representation, with x = [6. 9.6,6 C ] T , then the open loop 
equation can be expressed as 
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where 


x = Ax -I- Du; y = Cx; u = Qx + Rr 


The closed loop equation is expressed as 


x = Ax + B[Qx 4- Rr] = [A 4 bQ]x + K a br 
A -h hQ = 
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K a b = [0, 0, 0, (K a + K t )g] 

The state can be computed using Euler explicit integration 


x 


n+1 = At 


(A + bQ + —)x n + K a br n 


.n + 1 


= Cx 


n + 1 


The gain levels. K„ = 1.2 and K r = 0.06 s. were computed using conventional root locus methods. 
These gains were input into a linearized one degree-of-freedom dynamics routine to verify implementa- 
tion of the control law. The result of the linearized analysis is shown in Fig. 4. which shows a slightly 
divergent envelope for the canard-fixed case, and damped behavior for the closed-loop system. Note 
that a canard deflection limiter was applied to represent stalling of the airfoil. 
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Fig. 4: Linearized 2-D missile: body orientation and canard deflection time histories for both uncontrolled 
and controlled cases 


Nonlinear Coupled Simulation 


The control law developed above was also implemented into the code which coupled the RANS 
equations with rigid-body dynamics. From Fig. 3, the control law can be expressed as 


6c 


9 

s + f 


K a (9 c — 9) — K r 9 4" K t 0 c 


and integrated using Euler explicit integration 

*” +1 = f&t+ ' i { i Ka(6c ~ n ~ Kr ^ n + Kl9c \ + 

which can then be subjected to limit constraints. The result of the nonlinear solution shown in Fig. 5 
for both canard-fixed and controlled cases. After convergence to a steady-state solution, an accuracy- 
limited time step size of 46 ps was used, with a computational cost of five Cray Y-MP CPU hours. 
Grid communication was approximately one-third of the overall CPU time. Figure 5 shows similar 
behavior to the linearized cases, again with a modestly growing envelope for the canard-fixed case, and 
damped oscillations for the controlled case. This comparison provides a degree of validation of the 
implementation of the control law in the coupled code. 


2-D Resonating Cavity 

In order to establish confidence in the numerical method, a two-dimensional cavity computation 
was undertaken to demonstrate and validate self-induced cavity resonance. Details of this and three- 
dimensional cavity computations may be found elsewhere.' 

The test conditions, specified to match experiment, 22 were 

Moo = 0.9, Re L = 1.47 x 10 6 , L = 8 in. 

Poo — 0.40 kg/ iiv poc — 2.9 x 10 4 N jin 2 
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Fig. 5: 2-D missile: E 
for both uncontrolled ; 


Orientation 







and a cavity length t>y depth. L/D . of 2. Characteristic inflow and outflow conditions were specili* d 
a.li >ng with a step si of A/ = 1.97/z.s. Power spectra were computed using 8192 steps following the 
dissipation of tin.* initial transients. 

Inspection of the computed cavity pressure history, shown in Fig. Ga. confirms the idealized feedback 



Fig. 6: 2-d cavity: (a) pressure history and (b) power spectra comparison 

mechanism identified from Rossiter s experiments. Briefly, the cycle begins with the propagation of an 
acoustic wave from the aft wall of the cavity to the forward bulkhead. Wave reflection from the forward 
wall causes the shear layer to bow outwards, shedding vorticity. The deflected shear layer converts 
downstream and induces another cycle. This coupling of the acoustic and vortical fields is quantified 
by Rossiter's empirical model, given in Fig. 6b, which gives only feedback frequencies. 

In the frequency domain, comparison of Rossiter's data to present results indicate agreement in 
frequency at the peak magnitudes, as shown in Fig. 6b. Magnitudes are higher for the present case by 
about 2 dB. which can be explained from dimensionality arguments. The solution was also found to 
be insensitive to second-order dissipation levels within the range 0.3 to 0.5. Figure 6b also shows the 
resonant modes predicted by Rossiter's equation, showing that K = 0.56 gives better prediction of the 
higher modes. Finally, the vertical knife edge schlicren images of Fig. 7 show the qualitative agreement 
between computed and observed 23 radiation patterns. 
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Experiment, Karamcheti Present 


Fig. 7: Comparison of experimental and numerical schlieren images with knife edge vertical 

2-D Controlled Store Separation 

Using the control law developed for the 1-DOF simulation above, albeit with an arbitrarily chosen 
& c = 5°, a 3-DOF simulation of a store separating from a cavity was implemented. The separation 
began with the application of an 18000 N ejection force for 0.044 s , with corresponding acceleration 
of about 20g, during which the controller was off and no angular velocity was imparted to the store. 
The solution was initialized with the store fixed in carriage position, allowing damping of the starting 
numerical transients. Following convergence of the cavity acoustic envelope, the same accuracy-limited 
time step size of 46/i.s was used. The time step size was chosen such that the streamwise Courant 
number was about unity in the shear layer. This restriction is equivalent to allowing an acoustic wave' 
to propagate only one cell in a single step. Ten grids were used for this 1.3 x 10° point domain. The 
result, shown in Figs. 8 and 9, shows that the nose of controlled store remains pointed away from the 
parent body, while the canard-fixed store is pointed towards the parent 0.3 s after release. Since the 
controlled store is commanded to point away from the parent, the separation is faster for the coni ; - died 
case than for the canard-fixed store. Inspection of the normal force history shows a component of about 
50 Hz. corresponding to the second stage of Rossiter's formula. 9 


3-D Missile Stability Derivatives 

In order to determine the proper feedback gains, the stability derivatives of the missile must be 
computed. For this 3-DOF simulation this includes C ma , C rn ^ , and C mq . These parameters wore 
computed from four cases: using a nominal 0 = 8 C = 0. a constant pitch attitude 0 = 0.04 rad , 6 C = 0, 
a constant deflection angle 0 = 0, 8 C = 0.04 rad, and a constant pitch rate 0 = 9 rad/s . Canard motion 
was permitted by the yg in. nominal gap between the missile body and canard. This is a similar 
arrangement to the actual missile, albeit without the connector pin in the numerical model. The force 
and moment history was converged three orders of magnitude, approximately 2000 steps, on this 18 
grid solution containing approximately 1.5 million points. Figure 10 shows the geometry used for this 
portion of the study, along with coefficient of pressure, C p , contours. 

3-D Cavity Store Separation 

The geometry, shown in carriage position, can be seen in Fig. 11. The domain contains about 2.2 
million points distributed in 25 grids, with the missile grids being re-used from the previous stability 
derivative study. An example of the initialization of the cavity store separation problem is shown 
in Fig. 12. This store-fixed simulation will be run approximately five characteristic times, until the 
artificial starting transients have dissipated, after which ejection will begin. The ejection forces applied 
to the store will be such that the velocity at the end of the 8 in. piston stroke was 30 ft/s normal to 
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Fig. 8: 2-D controlled store separation: instantaneous Mach contours 
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Fig. 9: 2-D controlled store separation: time histories of a) center of gravity, b) force and moment 
coefficients, and c) body attitude and canard deflection angles 



Fig. 10: 3-D missile: C p contours 
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tho cavity with a pitch-down rate of 1 rad/s. These rates are the nominal conditions used during wind 
tunnel testing. 



Fig. 11: 3-D store separation: grids for missile in carriage 
Conclusions 

A pitch attitude control law was implemented in a coupled Navier-Stokes/rigid-body dynamics code. 
Comparison of nonlinear with linearized results showed reasonable comparison for the perturbed case 
with the controller active or off. Application of the control law to a two-dimensional, threc-degrec-of- 
freedom cavity store separation revealed improved trajectory characteristics. The generalized coding of 
aerodynamic effector kinematics in the coupled code will allow rapid implementation of existing control 
laws. Simulation of the coupled nonlinear aircraft trajectory can then be used to computationally 
prototype the control system. 


Note to the Reviewer 

Computation of a three-dimensional canard-fixed store separating from a cavity is in progress and 
will be compared to the experimental data of Dix and Dobson. 6 Comparison of the numerical and 
experimental results will provide an additional measure of validation. After completion of the uncon- 
trolled case, a simulation of the pitch-controlled case will commence. Determination of the feedback 
gains for this three-dimensional simulation will require the computation of the stability derivatives, 
which will also be accomplished via the Navier-Stokes equations. Assessment of the uncontrolled and 
controlled cases will show if improved cavity store separation characteristics can be achieved via canard 
effectors. 
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Fig. 12: 3-D store in carriage: instantaneous contours of pu on the symmetry plane and C p contours on 
the wetted surfaces 
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